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, We review here a novel circuit theory of superconductivity. The existed circuit theory of Andreev 

' reflection has been revised to account for decoherence between electrons and holes and twofold nature 

of the distribution function. The description of arbitrary connectors has been elaborated. In this 
way one can cope with the most of the factors that limited applicability of the old circuit theory. 
We give a simple example and discuss numerical implementation of the theory. 



43 



S 

I 



> 
in 
in 



oo 



I. INTRODUCTION 



Superconductivity is by virtue a mesoscopic phenomenon since the typical length scales involved exceed by far Fermi 
wavelenght and, frequently, the mean free path, jlj The adequate semiclassical theory of non-equlibrium superconduc- 
tivity has been initially elaborated for bulk systems. |^|| Subsequent contributions ||-[| that treated supeconducting 
constrictions and interfaces have provided the extension of the theory to heterogenious structures and the systems 



made artificially. Thereby the theoretical development has been essentially accomplished. The concise account for 

S this work one can find in |Q . We will refer to this patch of theoretical development as to " full theory" . 
There used to be a sharp contrast between the complexity of the full theory and relative simplicity of the final 
answers. It is a kind of disappointment to address, for instance, the problem of linear conductivity of a double tunnel 
barrier that separates a normal metal and a superconductor, to spend weeks and weeks calculating and to obtain that 
under very general assumptions G — Gn/^/2, Gn being conductivity in the normal state. 

The so-called " circuit theory of Andreev reflection" was born in an attempt to enhance the applicability of the full 
theory Initially it was a purely pedagogical project. I wanted to derive a primitive model that illustrates the essential 
features of the full theory. To my surprise, the model worked in a very efficient way. So I could not resist temptation 
and published a scientific article. || 

Since then, the circuit theory has been reviewed [p|Jl^1 and improved. C. W. J. Beenakker, D. Esteve, M. Devoret, N. 
Argaman reported extensions and reformulations, published as well as unpublished. Very elaborated formulation 
can be found in [Q . The author has extended the circuit theory to cover ballistic point contacts that enabled me to 
make calculations reported in 

In the present work we present a novel formulation of the circuit theory that incorporates most of the changes 
proposed. This theory enables energy- and voltage-dependent transport calculations in superconducting and normal 
structures. An important extension reported here is the concept of an "arbitrary connector". This allows for unified 
treatment of superconducting junctions of arbitrary nature. This extension is also relevant for the full theory. 

In fact, the resulting circuit theory looks very much like as a discrete version of the full theory. So that the full power 
■ of the full theory can now be incorporated in a circuit calculation. From the other hand, even a very sophisticated 
experimental layout can be presented with a few circuit theory elements. This essentially simplifies any practical 
I 1 calculation, numerical as well as analytical one. 

The outline of the paper is as follows. In Section || we remind the reader the approach of the old circuit theory 

o 
o 



and its limitations. We formulate the requirements to a circuit theory that is to overcome these limits. In Section III 
we provide a discrete version of the full theory equations in the diffusive limit. We introduce an important concept 
of the "leakage current" to describe decoherence between electrons and holes. We extend the concepts of a node and 
a connector in the next section. Section [v| is devoted to the arbitrary connector. We show there how a connector 
with known transmission eigenvalues shall be incorporated into the circuit theory. In section VI we accomplish the 



formula tion by giving the set of the rules of the new circuit theory and provide a simple example in Section VII 
Section VIII is devoted to the numerical algorithm that seem to suit most the structure of the theory. We discuss its 



practical implementations. 



II. CIRCUIT THEORY: OLD AND NEW 



Physics of electric conduction in a system that consists of superconductors and normal metals, with optional tunnel 
junctions, can be adequately treated with non-equilibrium superconductivity equations. Those are written for Keldysh 
Green's functions . (||| Those are 4x4 Nambu-Keldysh matrices matrices depending in a stationary case on space 
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coordinates and energy. They can be subdivided onto 2x2 Nambu matrices made up from usual and anomalous 
Green functions (we use 'check' for 4 x 4 and 'hat' for 2 x 2 matrices): 

C(t t>) (G R M G(x,x>) \ , A , ( G A (x,x') F A (x, x') \ 



G A (x,x / ) J ' 1 ' > ~ \-F +A (x,x') -G +A (x',x) 

and similar for G R and G. Semiclassical and diffusive appoximations allow to obtain equations for Green functions 
in coinciding points. We define G(x) — iv/irG(x, x), 

Here advanced and retarded functions determine the characteristics of energy spectrum of the system in a given point 
whereas K sets particle distribution over these energy states and thus directly related to the electric current and 
other physical quantities. It is assumed that the size of the system at least in a transport direction greatly exceeds 
Fermi wavelength and elastic mean free path, this makes diffusive approximation sensible. In this approximation, 
non-equilibrium superconductivity equations resemble a standard diffusion equation and read 

J^(v(x)G^G)-i[H,G]=0-, (3) 



dx a \ y ' dx° 

H = ea z + ia x Re(A(x)) + ia y lm(A(x)) . (4) 

Here A is superconducting pair potential, T> stands for diffusivity. A pseudounitary condition holds for the matrices 
G(e): G 2 = 1. We introduce here a density of a matrix current defined by 

j a (x)=a(x)G-^G, (5) 

a being specific conductivity in the normal state, so that the equation ^ can be presented as a conservation law for 
this current 

JLj<*(x)-ie 2 v[H,G]=0; (6) 

The matrix current density (||) can be related to the electric current density by means of 

jZ(x) = ^JdeTr[a z j K (e,x)] (7) 

If there are tunnel interfaces in the structure, G are in general different on both sides of the interface. Eq. [| shall 
be supplemented with the boundary conditions at the interface. These conditions have been derived in |6| and can be 
written in a comprehensive form as 

N ^ r{x) = 9M { G l{x) ^ 2 {x% (8) 

N a being a vector normal to the interface at a point x, g(x) being the conductance of the interface per unit area at 
the same point, 1,2 refer to different sides of the interface. These equations shall be also fulfilled by the boundary 
conditions "at infinity", otherwise their solution is not unique. Since we are talking about the transport, there must 
be at least two "infinities" that correspond to source and drain. In general, the structure can be connected to many 
bulk electrodes, either normal or superconducting, so it can have many "infinities". Such bulk electrode we will call 
" terminal" . In each terminal G shall assume its equilibrium value corresponding to the voltage at which the terminal 
is biased. 

With all these additions, Eq. [| provides a solid framework to calculate electric properties of superconducting 
structures. 

The circuit theory of Andreev conductance Q has grown from the fact that, under certain conditions, it is plausible 
to omit the second term in Eq. |^. So that the matrix current is conserved exactly. 

Now we spell limitations under which it is a true thing to do so. The first term is of the order oiV/L 2 , L being the 
system size, or more precisely, the size at which the resistance of the structure is being formed. In the normal metal 
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A = and the right term is of the order of e. Thus we need a sufficiently small system: L « ^V/e. It implies 
that: a) the temperature is low enough: T < V/L 2 , T < A. b) the voltage is low enough: eV < V/L 2 . V < A. The 
last limitation below is given by the fact that we use stationary equation. If there had been several superconducting 
terminals in the structure biased by different voltages, it would have given rise to non-stationary Josephson-like effect 
that would make the Green functions to depend on two energies. So that, c) all sureconducting terminals are at the 
same voltage. Let us set this voltage to zero. 

The key idea of any circuit theory is to get from continuous conservation law of the type (Q) to its discrete version. 
The structure is presented as a set of nodes those are connected pairwise by means of connectors. The conservation 
law is presented as a set of second Kirchhoff rules for each node expect terminal ones, 

= ^4, (9) 

k 

where the summation goes over all the nodes k connected to the node i. The current Iik — —Im depends on the states 
of the nodes i, k and on the connector. Therefore, the Eqs. || determine the state of the nodes at a given state of the 
terminals. In a normal terminal biased at voltage V G is given by 



2 r> /tanh^ . 

R =- A = °- K -\ -tanh— I' 



In a superconducting terminal, 
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hA = h*R = H; K = (R - A) tanh (11) 



where h = y/ (e — iS) 2 — | A| 2 . In a common electric circuit theory the state of the node is characterized by its 
electrostatic potential. The current is a scalar given by Iik — gik{4>i ~ 4>k)- In the circuit theory || the state of the 
node is in principle characterized by the full G(e) and the current retains matrix structure. However the symmetries of 
G specific for e ~ allow for significant simplifications. Energy dependence of A, R can be disregarded and K can be 
integrated over the energy. The Kirchhoff equations can be separated on two parts. First part are equations for vector 
"spectral currents" to determine "spectral vectors" A, R in each node. Second part are equations for usual currents 
to determine non-equilibrium chemical potentials in each node and, finally, Andreev conductance. The resistance of 
the tunnel connectors is remormalized, this renormalization is given by scalar product of two "spectral vectors" on 
the sides of the connector. 

The old circuit theory is good to obtain simple answers. If the analysis of experimental data in real structures is 
meant, the conditions a) and b) appear to be very restrictive. Beside the fact that the temperature can not be made 
arbitrary low, the voltage and temperature dependence of the transport can not be accessed in the framework of the 
theory, j§ 

We will not explain here the rules and details of the old circuit theory: those can be found in f§||. Our goal is 
to formulate a circuit theory that is free form e ~ limitation. We can now spell the requirements to the structure 
of this theory. First, in order to comply with the full theory, the state of the node shall be characterized by a full 
G matrix. Second, since energy dependence of the Green functions has to be taken into account, the equations (|^) 
shall be written and solved separately in each energy slice. By doing so we disregard possible inelastic scattering, 
which is usually not important. Physical values like electric current are then obtained by integration over e. The 
third requirement stems from the recent technological achievements. It is possible now to incorporate (ballistic) point 
contacts into mesoscopic superconducting structures. So we require that the theory shall describe such contacts and, 
in general, any arbitrary constrictions and connections. 

Below we derive the rules of this theory and discuss their implementations. 



III. DISCRETIZATION OF DUFFUSIVE CONDUCTOR 

Let us begin with formulation of the discrete version of Eq. |]. Instead of continuous space x a we take a connected 
discrete set of xf such that G(xf) in neighboring points of the set are close to each other. We associate a resistor 
with each connection in the set in such a way that it simulates continuous resistivity of the system. Let us show how 
to choose such resistors for square lattice of xf with periods a" 2 ( a^af = a^a^ — a 2 ,afa,2 = 0) that approximates 
homogeneous film with sheet conductivity a. Let us expand G in the vicinity of the node i: G(x) = G(xi)-\-'E a (x a —xf), 
S"a << G(x,i). From Eq. p] we obtain the continuous matrix current density 
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r = aG( Xi )E a (12) 

In the network, the current density is given by 

r=a- 2 J2^K-xf) (13) 

k 

k being nodes neighboring i. We see that if we choose 

Iik = 9 -f[GuG k ], (14) 

gik = era we reproduce Eq. |l2]. To prove this we make use of the fact that [E, a ,Gi}+ = 0, the latter follows from 
G 2 = 1. Now we can rewrite Eq. || as 

= Iik + Ileakage', (15) 

k 

where we introduce "leakage" current 

Ileakage = -ie 2 vV l [G i , H], (16) 

Vi being volume associated with the node i. The term "leakage" comes from the fact that the Eq. has a formal 
similarity with the equations of the electric circuit theory that take into account the charge leakage to the ground 
from each node. There is no net leakage of the electric current due to ([if]) since the corresponding term appears to 
be zero. However, the matrix current (|l^) describes two processes that may be viewed as a sort of leakage. Namely, 
the terms proportional to e describe decoherence between electrons and holes, that is, the fact that the electrons and 
holes at the same energy difference e from Fermi surface have slightly mismatching wavevectors. So we have "leakage 
of coherence" . The terms proportional to A are responsible for the conversion between quasiparticles and cooper 
pairs that form the superconducting condensate. This is a leakage of quasiparticles. 

We note that, by virtue of Eq. |], the Kirchhoff rules (|l4|), ( |i"5| ) are also suitable for systems containing tunnel 
junctions. Just some connections of the resistive network shall be replaced with (specific) conductances of the tunnel 
barriers. 

With increasing fineness of the node set we obviously converge to continuous limit and obtain more and more 
accurate agreement with the results of the full theory. A practical question is how fine should this mesh be in a realistic 
calculation. A short answer is that the nodes shall be close than than the coherence length £ = ^JV max(A, e). In 
practice, the mesh describing a diffusive wire should contain about a dozen nodes even at £ that exceeds the system 
size. This drawback is compensated in the limit of short £, since in this case the interesting behaviour occurs merely 
at the edge of the system and it makes no sense to keep fine mesh through the whole sample. It turns out that the 
mesh of 40 nodes is sufficient to describe transport properties of the diffusive wire this the accuracy at least 10 -4 in 
the whole energy range. 



IV. MORE ABOUT NODES AND CONNECTORS 

As we saw above, the essential requirements to a node are i. its state is characterized by G, so that Green functions 
must be isotropic on the scale on the node size, ii. the size of the node is less than coherence length. 

This allows us to treat as nodes parts of the system where the transport is not or is not entirely diffusive. Let us 
consider, for instance, the well-known model of ballistic cavity, jlij Although the transport in the cavity is the ballistic 
one, one can actually regard the cavity as a node connected to reservoirs by means of ballistic point contacts. It is the 
chaotic character of ballistic transport that makes Green function isotropic. The condition of good isotropization is 
fulfilled provided the diameter of the cavity greatly exceeds the typical size of the point contacts, or, in other terms, 
the conductance of the contacts must be much less than the Sharvin estimation of the cavity conductance. The cavity 
must be smaller than the coherence length in the ballistic limit. The leakage current from the cavity is still given by 
Eq. it does not know which mechanism has provided the isotropization of the Green function, the cavity could 
be diffusive as well. 

This makes it relevant to investigate such connectors between the nodes as ballistic point contacts. Previous 
experience with superconducting constrictions shows a significant difference between tunnel junctions, diffusive and 
ballistic constrictions. This difference stems from the fact that the constrictions of the same conductance have different 
distribution of transmission eigenvalues. 

This implies that the expression ( |l4| ) is not valid for the arbitrary connector. Below we derive the expression for 
the matrix current in terms of transmission eigenvalues. 
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V. ARBITRARY CONNECTOR 



The derivation presented below makes substantial use of the approach || and the general scattering theory outlined 
in fTs] ]. The goal is to express the matrix current in the junction between two nodes in terms of isotropic Green 
functions in these nodes and transmission eigenvalues that characterize the scattering in the junction. We subdivide 
the system onto five zones: two isotropization zones on both sides of the junction, two ballistic zones and scattering 
zone. We assume that the size of all zones in transport direction is much shorter than the coherence length. Under this 
condition, the matrix current conserves and can be easily evaluated in ballistic zone where Green function equations 
are simple. The same condition allows us to disregard energy dependence of the scattering amplitudes. 

Following to |l5| , we discretize the transverse motion of electrons. We introduce a number of transport channels 
labeled by n and let z to be along the transport direction. The scattering takes place near z = 0. Let us first consider 
the left side of the junction, z < 0. Similar to the approach of the ref. j|] we present the exact Green function 
G n m{z, z') in the following form 

G nm (z,z')= ^ exp(i<7Pn-i<r'Pm)Gnm{z,z'), (17) 

a,cr'=±l 

p n being Fermi momentum in the n-th channel. The functions G aa nm are varying smoothly at the scale of l/p n and 
obey the following semiclassical equations: 

(<Jv n §- z +H)G™ m =0 (18) 
^m(^«^7-#) = 0. (19) 

Ti. being the effective Hamiltonian that comprises H and the self energy part that describes isotropization scattering. 
These functions are discontinuous at z = z', their values for z > z' and z < z' being matched with the aid of the 
following condition [|| 

GZ n (z + 0,z)- GZa (z-0,z) = iia8 nm 5 aa , (20) 

It is convenient to work with matrices defined in an extended basis with comprises Nambu-Keldysh indices, channels, 
and a= direction of the mode. We will use tilde accent for matrices in this basis, for example, A We will also denote 
with A the matrices that are diagonal in mode indices and with A those diagonal in Nambu-Keldysh indices. To 
characterize the Green functions we introduce the matrix g such that 

2i6°£(z,z') =^(z,z f )/ y /w^ + (r8 mr >siffi(z-2l , )/\v n \ (21) 

The matrix g is continuous at z — z'. The factors v n are chosen in such a way that g is composed of the scattering 
wave functions that are normalized per unit flux. This makes it easy to apply the general scattering formalism of Ref. 
p5[ . Following we introduce a transfer matrix that relates wave functions on the right side (ip^) and those on 
the left side (V^ 1 ^) of the scatterer, 

^11 = E M m,*;n,*^ { nl> (22) 
n,(7' 

Since g transforms as a product of wave functions, its values at z, z' > and z, z 1 < are related by 

g 2 = M f ~ 9l M. (23) 

The transfer matrix M obeys flux conservation , 

ME z M f = E z , (24) 

where m n = <7& m nbaa'- It follows from ( ]2~I| ) and ( p3| ) that the quantity we are after, the matrix current, can be 
expressed in terms of g on either side: 

/ = —Tr m ,„ [t'gt] = -Tw [t z g 2 ] . (25) 
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To find the matrix current, we have to evaluate gin)- To do so, we shall consider the behavior of g in the isotropization 
zone on both sides of the scatterer. Let us first consider the left side of the scatterer and let us assume the simplest 
model of the isotropization: scattering on the point defects. Then in the isotropization zone we may approximate 



Iri = Gi/2Ti mp . Since H does not depend on z the equations 19 can be readily solved in the most general form and 



we can express g(z, z') in the isotropic zone in terms of its value g\ in the ballistic zone, 

G^(z, z') = P(z)(~ gi + sign(z - z')^)P(-z'); (26) 
where P is diagonal over channel index, 

P{z) = ^(expt— -— )(1 - + exp(-— + E'fix)) (27) 

To derive that we use that G\ = 1. To assure that Gt^m does not grow with decreasing z, we shall require 

(£*+<2i)(£*-0i)=O; (28) 

(£ z +.gi)(£*-Gi) =0. (29) 
Similar consideration for the right side of the scatterer yields 

(t z -G 2 )(t z + ~g 2 )=0; (30) 

(E z -g 2 )(£* + G 2 ) =0. (31) 

We see that the above conditions do not contain any information about how the isotropization actually occur. This 
suggests that those are of universal nature and do not depend on isotropization mechanism. 

Conditions (p8|)-(|3l|) together with Eq. |23| co mpletely determine <?i( 2 ), Green function in the ballistic region. To 
find an explicit expression, we multiply Eq. by M' from the left and by M from the right. Taking into account 
Eqs. |2^ and we find that 

Q + GiE 2 -(QE z + Gi).g 2 = 0. (32) 

Following [jl5|| , we have introduced here the matrix Q = M^M, Q _1 — Y^ Z QY7- . Now we multiply Eq. ^ by Q from 
the left and add it to Eq. [52]. We derive g 2 from the resulting equation, 

.92= ^ * (2Q + (G 1 -QG 2 )S Z )). (33) 



After some algebra it is possible to check that expression (^3|) satisfies the conditions (|28|)-(|3l|) provided G\ = G\ — 1, 
as it should be. Besides, <?i( 2 ) = !■ The equation |33| presents an important step in the derivation. 

To evaluate the matrix current we find g 2 in the basis composed of eigenvectors of Q and Nambu-Keldysh indices. 
In this basis g 2 takes a blockdiagonal form. It appears Jl5| that for each eigenvector c„ with eigenvalue q n > 1 the 
vector T, z c n is also an eigenvector of Q with eigenvalue q^ 1 . These two eigenvectors form a block. Within the block 

Qn V *z _ ( 1 



This allows us to find g 2 within the block: 



{Gx-qGi) 



■92 = i 9 Z (a F C \ qM "2^__ ■ (35) 



Now we can directly apply Eq. 25 to evaluate the matrix current. Finally we obtain 



— IT G 2 Gi - GiG 2 

n "4 + T„(G 2 G 1 +G 1 G 2 -2) 



I = — 2T» ; - ^1 V't (36) 
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where we have switched to more conventional representation of q n in terms of transmissions T n , eigenvalues of the 
transmission matrix square. 

Taking its practical use apart, the relation ( |36| ) is probably the most concise way to express the accumulated 
knowledge about superconducting constrictions. We note that the electric current in the constriction is given by the 
Keldysh part of (EF 



hi = — Tr 

4e 



(37) 



If we set one of G to the values of a normal terminal ( |lC)| ) and another one to the values of superconducting terminal 
( |TT| ) we obtain i. Beenakker's relation for Andreev conductance in the limit e — > 0, ii. a simple generalization 
of Blonder-Tinkham-Klapwijk formula Q at e of the order of A. Setting both G to superconducting terminals of 
different phases give all known relations for the Josephson current. 

Keldysh components of (^6|) can be separated. For retarded one we have 

JR = i_Y — ^— r[iMl] (38) 

Tr ^2 + T n Ri-R 2 -1 y 1J K ' 

and similar for A. This suggest that we can characterize the scatterer with a single scalar function Z(x), 

n v ' 

For any constriction Z(l) = Gat/2. 

The expression for Keldysh component looks rather cumbersome: 

I K = Z(Ri • R 2 ) J R 1 i^ 2 + Z( R '' R ' } ~ Z A (Al A A2) hk 2 A 2 + R 1 K 2 A 1 A 2 - PkkxkiA^) - Z(A 1 ■ A 2 )K 2 A 1 - (1 <-> 2) 



2(Ri ■ Ri - Ai ■ A 2 



(40) 



It is conventional to present K in the form Rf — fA, where / = f L l + f T a z presents two-component distribution 
function in the superconductor. We note that the expression ( |40| ) has in fact only two independent components, that 
corresponds to two components of the distribution function. One can choose these components to be Iq — TiI K , 
I z = Tra z I K . Then one can express these components in terms of f T,L on both sides of the junction. Later one can 
keep only these components in the Kirchhoff rules ([l5|). However, the corresponding expressions are so cumbersome 
that we do not dare to recommend this way: it is simpler to work with all four components of I K , 

Albeit there is an important case where the presentation mentioned helps a lot. If there is no mixing of normal and 
superconducting current whithin the system, equations for Io and I z separate, with I z depending on f T only, 

2I z (e)=G(e)(f 2 T (e)~f?(e)); (41) 

where G(e) reads 

G(e) = 2Fie (Z(00 + 2ReZ(£)Re M + ^"J^ (Im| 2 + l£| 2 + 2Re/iRe£ - (1 - fll )(l - s 2 )) (42) 

where 

A1 = A+-A 2 ; C = A+-A+; s 1)2 = A+ 2 • Ai, 2 . (43) 

It is important that G(e) can be viewed as energy-dependent effective conductance of the junction. This follows from 
the fact that the Kirchhoff rules that determine f T at any given energy can be viewed as electric circuit theory rules 
with energy-dependent conductances. 
For tunnel junctions, Eq. ^2| reduces to 

G(e) = ^Tr [(ii + A+)(i 2 + A+)] (44) 

The "no mixing" condition is always satisfied if there is only only superconducting terminal in the system. It can 
be satisfied if the normal and superconducting current are geometrically separated like in the experiment. [fL6| 
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VI. RULES 



We summarize the above derivation by giving the rules of the novel circuit theory. 

i. Define the circuit. That includes a proper choice of the node mesh in the diffusive parts of the system, determi- 
nation of the type of the connectors in terms of Z(x). 



ii. Write down Kirchhoff rules (15) with leakage current. Information about connectors shall be used at this point. 

iii. Setting terminals to given values (|l^Jll|) , find from the Kirchhoff rules G in each node for each energy. 

iv. Find currents in the circuit with the aid of Eqs. |io| , p^ ,[7| 

These rules are to apply in the most general case. If "no mixing" condition is satisfied, Kirchhoff rules have to be 
solved for R, A only. The energy-dependent conductances of all connectors are then given by Eq. |42. 

In the rest of this section we give a simple relation that allows practical work with Eqs. [l5] and occasionally gives 
analytical results. The point is that by virtue of (pq) Iik can be always presented in the form of Iik — [Gi,Lik], L 



being a matrix. So that the second Kirchhoff rule for the node i takes the commutator form 

=[&,&]. (45) 
In the simple case where all connectors are of tunnel nature expression for Q\ reads 

Gi = Y d 9ikG k + evV i H (46) 

k 

and does not depend on Gi . In general, it depends on G; by virtue of Eq. |3(| In any case the Eq. allows to express 
G in terms of Q\. We explicitly write Keldysh components of these matrices: 

Moi>=(o5> <«> 

For diagonal blocks we find that 

l = i/a; R = n/r; (48) 
where a = y/Tr(A 2 )/2, r = a* . Substituting this into ( f45| ) we obtain that 

k = {sL + r)- l {K-RKA). (49) 



VII. SIMPLE EXAMPLE 



To give a simple application of the circuit theory rules we consider a system of two tunnel junctions of the same 
resistance Rt with normal metal in between that separate a superconducting electrode and a normal metal biased 
at voltage V. If we assume that the correlation length is shorter than the size of the intermediate normal metal and 
its resistance is negligible in comparison with the resistances of the tunnel junctions we can regard it as a node. This 
node is connected to two terminals. With the aid of (pi]), (10), (p^) , (f48|) we find the Green function in the node. To 
simplify the formulas, we write it explicitly only for e <C A: 



1 + ier 1 



+ ier) 2 + 1 V 1 -1-ier 



(50) 



where r = e 2 vV no d e RT is a typical escape time from the node. Since there is no superconducting current in the 
system, we do not have to write K explicitly. Instead, we evaluate the energy-dependent resistances for each junction 
with the aid of Eq. || 

Rt = Rt r- ;Jfr = R T -y - (51) 

Re (1/^(1 -Her) 2 + lj Rc (Jl + ier)/y/(l + ier) 2 + lj 

The differential conductance G(V) of the system is then given by l/(Ri + R2) at e = eV if T = 0. Fig.l presents 
plots of R T /Ri,2,G(V). 
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VIII. NUMERICAL IMPLEMENTATION 



We have developed a simple numerical algorithm that suits the structure of the circuit theory proposed. It is based 
on successive iterations of Eqs. 48 4^ at a given energy. Initially, an approximation for G in each node is stored in 
memory Next step is calculation of Q for each node. Here we make use of information about connections in the circuit. 
From Eqs.^,^ we obtain then next approximation for G,. This algorithm is not expected to converge very fast. In 
fact, the calculation time is proportional to N^ odes . However, the algorithm does not make use of any complicated 
parametrization of G, that speeds up the calculation. Also, Green functions usually have to be found at many close 
energies and are calculated consequently. In this case, one can use the final result for G at a certain energy as a good 
initial approximation for the next energy. All this makes the algorithm fast enough to suit the practical purposes. 
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